✅ Key points
- Opening line: starts with three backticks +
{r ...}- Contains the chunk name (
chunk_name) and options (echo,eval, etc.)
- Contains the chunk name (
- Code: written below the opening line
- Closing line: three backticks alone
(
\```)- Nothing else should be on this line
- Nothing else should be on this line
- Optional: you can add comments inside the chunk (
#) to describe what the code does
| Option | Purpose |
|---|---|
echo=TRUE |
Shows the R code |
eval=TRUE |
Runs the code (otherwise nothing is executed) |
message=FALSE |
Hides package startup messages |
warning=FALSE |
Hides warnings in the HTML |
complete_2023_2024 <- read.csv("input/2023_2024_all_moth_counts.csv")
#clean_complete - full summer complete moth counts available for 2023 and 2024
library(tidyverse)
library(lme4)
library(performance)
library(summarytools)
library(ggeffects)
library(marginaleffects)
library(brms)
library(car)
library(viridis)
library(broom)
library(knitr)
library(dplyr)
library(kableExtra)
library(gt)
library(modelsummary)
library(sjPlot)
# Data Cleaning -----------------------------------------------------------
##To explore the distribution of your variables and count data like clean_complete
# quick visualizations
dfSummary(complete_2023_2024)## Data Frame Summary
## complete_2023_2024
## Dimensions: 255 x 20
## Duplicates: 0
##
## ------------------------------------------------------------------------------------------------------------------------
## No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
## ---- ------------------- ----------------------------- --------------------- --------------------- ---------- ----------
## 1 patch_name 1. Oka 44 (17.3%) III 255 0
## [character] 2. Orford 39 (15.3%) III (100.0%) (0.0%)
## 3. Rigaud 30 (11.8%) II
## 4. Mont_Saint_Hilaire 29 (11.4%) II
## 5. Montebello 22 ( 8.6%) I
## 6. Notre_Dame_de_Bonsecours 18 ( 7.1%) I
## 7. Mont_Royal 13 ( 5.1%) I
## 8. Mont_Saint_Bruno 12 ( 4.7%)
## 9. Papineauville 12 ( 4.7%)
## 10. Brownsburg 6 ( 2.4%)
## [ 5 others ] 30 (11.8%) II
##
## 2 Years Min : 2023 2023 : 102 (40.0%) IIIIIIII 255 0
## [integer] Mean : 2023.6 2024 : 153 (60.0%) IIIIIIIIIIII (100.0%) (0.0%)
## Max : 2024
##
## 3 Year Min : 0 0 : 102 (40.0%) IIIIIIII 255 0
## [integer] Mean : 0.6 1 : 153 (60.0%) IIIIIIIIIIII (100.0%) (0.0%)
## Max : 1
##
## 4 stand_type 1. MOM 10 ( 3.9%) 255 0
## [character] 2. Oak 73 (28.6%) IIIII (100.0%) (0.0%)
## 3. Oak/Other 33 (12.9%) II
## 4. Oak/Pine 41 (16.1%) III
## 5. Other 30 (11.8%) II
## 6. Pine 39 (15.3%) III
## 7. Pine/Oak 29 (11.4%) II
##
## 5 landscape_type 1. agricultural 139 (54.5%) IIIIIIIIII 255 0
## [character] 2. forest· 85 (33.3%) IIIIII (100.0%) (0.0%)
## 3. urban 31 (12.2%) II
##
## 6 stand_category 1. C 35 (17.9%) III 196 59
## [character] 2. I 29 (14.8%) II (76.9%) (23.1%)
## 3. E 24 (12.2%) II
## 4. H 18 ( 9.2%) I
## 5. F 17 ( 8.7%) I
## 6. L 15 ( 7.7%) I
## 7. A 10 ( 5.1%) I
## 8. D 10 ( 5.1%) I
## 9. K 10 ( 5.1%) I
## 10. MOM 10 ( 5.1%) I
## [ 4 others ] 18 ( 9.2%) I
##
## 7 Percent_Oak Mean (sd) : 0.3 (0.2) 0.00 : 48 (18.8%) III 255 0
## [numeric] min < med < max: 0.10 : 22 ( 8.6%) I (100.0%) (0.0%)
## 0 < 0.3 < 0.8 0.20 : 10 ( 3.9%)
## IQR (CV) : 0.4 (0.7) 0.30 : 58 (22.7%) IIII
## 0.40 : 50 (19.6%) III
## 0.50 : 18 ( 7.1%) I
## 0.60 : 18 ( 7.1%) I
## 0.70 : 19 ( 7.5%) I
## 0.80 : 12 ( 4.7%)
##
## 8 Percent_Pine Mean (sd) : 0.2 (0.2) 0.00 : 117 (45.9%) IIIIIIIII 255 0
## [numeric] min < med < max: 0.10 : 29 (11.4%) II (100.0%) (0.0%)
## 0 < 0.1 < 0.8 0.20 : 27 (10.6%) II
## IQR (CV) : 0.4 (1.2) 0.30 : 12 ( 4.7%)
## 0.40 : 21 ( 8.2%) I
## 0.50 : 19 ( 7.5%) I
## 0.60 : 9 ( 3.5%)
## 0.70 : 15 ( 5.9%) I
## 0.80 : 6 ( 2.4%)
##
## 9 Percent_Maple Mean (sd) : 0.3 (0.2) 10 distinct values II 246 9
## [numeric] min < med < max: III (96.5%) (3.5%)
## 0 < 0.2 < 0.9 IIII
## IQR (CV) : 0.3 (0.6) IIII
## IIII
##
## 10 Forest_Type 1. oak 65 (25.5%) IIIII 255 0
## [character] 2. oak_maple 48 (18.8%) III (100.0%) (0.0%)
## 3. pine 45 (17.6%) III
## 4. oak_pine 21 ( 8.2%) I
## 5. maple_oak 16 ( 6.3%) I
## 6. pine_oak 14 ( 5.5%) I
## 7. maple 12 ( 4.7%)
## 8. maple_birch 7 ( 2.7%)
## 9. maple_hardwood? 5 ( 2.0%)
## 10. pine_hemlock 5 ( 2.0%)
## [ 9 others ] 17 ( 6.7%) I
##
## 11 stand_area_ha Mean (sd) : 8.4 (3.5) 67 distinct values . : 255 0
## [numeric] min < med < max: : : : (100.0%) (0.0%)
## 3.5 < 8 < 25.9 : : :
## IQR (CV) : 3.2 (0.4) : : : :
## : : : : : .
##
## 12 forest_area_km2 Mean (sd) : 210.2 (436.8) 13 distinct values : 255 0
## [numeric] min < med < max: : (100.0%) (0.0%)
## 1.6 < 14 < 1282 :
## IQR (CV) : 95.7 (2.1) :
## : :
##
## 13 trap_name 1. BC Co-Dom1 1 ( 0.4%) 255 0
## [character] 2. BC Co-Dom2 1 ( 0.4%) (100.0%) (0.0%)
## 3. BC Dom1 1 ( 0.4%)
## 4. BC Dom2 1 ( 0.4%)
## 5. BC Low 1 1 ( 0.4%)
## 6. BC Low 2 1 ( 0.4%)
## 7. Bolt Co-Dom1 1 ( 0.4%)
## 8. Bolt Co-Dom2 1 ( 0.4%)
## 9. Bolt Dom1 1 ( 0.4%)
## 10. Bolt Dom2 1 ( 0.4%)
## [ 245 others ] 245 (96.1%) IIIIIIIIIIIIIIIIIII
##
## 14 longitude Mean (sd) : -73.7 (1) 255 distinct values : 255 0
## [numeric] min < med < max: : (100.0%) (0.0%)
## -75 < -74 < -72 : : : :
## IQR (CV) : 1.2 (0) : : : :
## : : : : : : .
##
## 15 total_moth_count Mean (sd) : 57.6 (51.4) 113 distinct values : 252 3
## [integer] min < med < max: : (98.8%) (1.2%)
## 0 < 47.5 < 378 : .
## IQR (CV) : 55.5 (0.9) : :
## : : :
##
## 16 complete Mean (sd) : 72.8 (60.9) 92 distinct values : 142 113
## [integer] min < med < max: : : (55.7%) (44.3%)
## 3 < 61 < 378 : :
## IQR (CV) : 68.8 (0.8) : : :
## : : : . .
##
## 17 clean_complete Mean (sd) : 62.3 (55.4) 140 distinct values : 200 55
## [numeric] min < med < max: : . (78.4%) (21.6%)
## 0 < 49.8 < 378 : :
## IQR (CV) : 59 (0.9) : :
## : : :
##
## 18 censored Min : 0 0 : 200 (78.4%) IIIIIIIIIIIIIII 255 0
## [integer] Mean : 0.2 1 : 55 (21.6%) IIII (100.0%) (0.0%)
## Max : 1
##
## 19 reason_incomplete 1. (Empty string) 84 (60.4%) IIIIIIIIIIII 139 116
## [character] 2. Missing in field 2 ( 1.4%) (54.5%) (45.5%)
## 3. Muck 44 (31.7%) IIIIII
## 4. Trap damaged· 9 ( 6.5%) I
##
## 20 X All NA's 0 255
## [logical] (0.0%) (100.0%)
## ------------------------------------------------------------------------------------------------------------------------str(complete_2023_2024)## 'data.frame': 255 obs. of 20 variables:
## $ patch_name : chr "Mont_Royal" "Mont_Royal" "Mont_Royal" "Mont_Royal" ...
## $ Years : int 2024 2024 2024 2024 2024 2024 2024 2024 2024 2024 ...
## $ Year : int 1 1 1 1 1 1 1 1 1 1 ...
## $ stand_type : chr "Oak" "Oak" "Oak" "Oak" ...
## $ landscape_type : chr "urban" "urban" "urban" "urban" ...
## $ stand_category : chr "E" "A" "C" "C" ...
## $ Percent_Oak : num 0.4 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.4 0.4 ...
## $ Percent_Pine : num 0 0 0 0 0 0 0 0 0.1 0.1 ...
## $ Percent_Maple : num 0.3 0.1 0.3 0.3 0.3 0.3 0.3 0.3 0 0 ...
## $ Forest_Type : chr "oak_maple" "oak" "oak" "oak" ...
## $ stand_area_ha : num 14 6.4 8.6 8.6 8.6 11.7 11.7 11.7 8.9 8.9 ...
## $ forest_area_km2 : num 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 ...
## $ trap_name : chr "MRMOM" "MROak2.1" "MROak1.1" "MROak1.2" ...
## $ longitude : num -73.6 -73.6 -73.6 -73.6 -73.6 ...
## $ total_moth_count : int 93 15 7 20 6 96 52 NA 60 20 ...
## $ complete : int 93 15 7 20 6 96 52 NA 60 20 ...
## $ clean_complete : num 93 15 7 20 6 96 52 NA 60 20 ...
## $ censored : int 0 0 0 0 0 0 0 1 0 0 ...
## $ reason_incomplete: chr NA NA NA NA ...
## $ X : logi NA NA NA NA NA NA ...## change 'Year' from int to chr
complete_2023_2024$Year <- as.character(complete_2023_2024$Year)
str(complete_2023_2024)## 'data.frame': 255 obs. of 20 variables:
## $ patch_name : chr "Mont_Royal" "Mont_Royal" "Mont_Royal" "Mont_Royal" ...
## $ Years : int 2024 2024 2024 2024 2024 2024 2024 2024 2024 2024 ...
## $ Year : chr "1" "1" "1" "1" ...
## $ stand_type : chr "Oak" "Oak" "Oak" "Oak" ...
## $ landscape_type : chr "urban" "urban" "urban" "urban" ...
## $ stand_category : chr "E" "A" "C" "C" ...
## $ Percent_Oak : num 0.4 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.4 0.4 ...
## $ Percent_Pine : num 0 0 0 0 0 0 0 0 0.1 0.1 ...
## $ Percent_Maple : num 0.3 0.1 0.3 0.3 0.3 0.3 0.3 0.3 0 0 ...
## $ Forest_Type : chr "oak_maple" "oak" "oak" "oak" ...
## $ stand_area_ha : num 14 6.4 8.6 8.6 8.6 11.7 11.7 11.7 8.9 8.9 ...
## $ forest_area_km2 : num 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 4.89 ...
## $ trap_name : chr "MRMOM" "MROak2.1" "MROak1.1" "MROak1.2" ...
## $ longitude : num -73.6 -73.6 -73.6 -73.6 -73.6 ...
## $ total_moth_count : int 93 15 7 20 6 96 52 NA 60 20 ...
## $ complete : int 93 15 7 20 6 96 52 NA 60 20 ...
## $ clean_complete : num 93 15 7 20 6 96 52 NA 60 20 ...
## $ censored : int 0 0 0 0 0 0 0 1 0 0 ...
## $ reason_incomplete: chr NA NA NA NA ...
## $ X : logi NA NA NA NA NA NA ...# remove space in trap_name
# replace '-' with '_' in trap_name
complete_2023_2024 <- complete_2023_2024 %>%
mutate(trap_name = str_replace(trap_name, " ", ""))
complete_2023_2024 <- complete_2023_2024 %>%
mutate(trap_name = str_replace(trap_name, "-", "_"))
#Remove MOM traps from either "stand type" or "stand category"
stand_type_filtered <- complete_2023_2024 %>%
filter(stand_type != "MOM")
stand_category_filtered <- complete_2023_2024 %>%
filter(stand_category != "MOM")
# looking for mistakes
unique(complete_2023_2024$stand_type)## [1] "Oak" "Oak/Pine" "Pine" "MOM" "Pine/Oak" "Oak/Other"
## [7] "Other"unique(complete_2023_2024$patch_name)## [1] "Mont_Royal" "Mont_Saint_Hilaire"
## [3] "Montebello" "Notre_Dame_de_Bonsecours"
## [5] "Oka" "Orford"
## [7] "Papineauville" "Rigaud"
## [9] "Brownsburg" "Hatley"
## [11] "Kenauk" "Mont_Gauvin/Glen "
## [13] "Mont_Saint_Bruno" "Parc_Michel_Chartrand"
## [15] "Parc_Pointe_aux_Prairies"unique(complete_2023_2024$Year)## [1] "1" "0"unique(complete_2023_2024$landscape_type)## [1] "urban" "agricultural" "forest "unique(complete_2023_2024$stand_category)## [1] "E" "A" "C" "L" "MOM" "D" "H" "I" "B" "G" "K" "M"
## [13] "F" "J" NAunique(stand_type_filtered$stand_type)## [1] "Oak" "Oak/Pine" "Pine" "Pine/Oak" "Oak/Other" "Other"unique(stand_category_filtered$stand_category)## [1] "E" "A" "C" "L" "D" "H" "I" "B" "G" "K" "M" "F" "J"# Data Summaries ----------------------------------------------------------
##separate Stand Type column so that we have a Stand ID for each stand in each patch
stand_ID_filtered <- stand_type_filtered %>%
separate(trap_name, into = c("stand_ID", "trap_number"), remove = FALSE, sep = "\\." ) %>%
glimpse()## Rows: 245
## Columns: 22
## $ patch_name <chr> "Mont_Royal", "Mont_Royal", "Mont_Royal", "Mont_Roya…
## $ Years <int> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024…
## $ Year <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1…
## $ stand_type <chr> "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oa…
## $ landscape_type <chr> "urban", "urban", "urban", "urban", "urban", "urban"…
## $ stand_category <chr> "E", "A", "C", "C", "C", "E", "E", "E", "E", "E", "E…
## $ Percent_Oak <dbl> 0.4, 0.8, 0.6, 0.6, 0.6, 0.4, 0.4, 0.4, 0.4, 0.4, 0.…
## $ Percent_Pine <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.1, 0.…
## $ Percent_Maple <dbl> 0.3, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.0, 0.0, 0.…
## $ Forest_Type <chr> "oak_maple", "oak", "oak", "oak", "oak", "oak_maple"…
## $ stand_area_ha <dbl> 14.0, 6.4, 8.6, 8.6, 8.6, 11.7, 11.7, 11.7, 8.9, 8.9…
## $ forest_area_km2 <dbl> 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89…
## $ trap_name <chr> "MRMOM", "MROak2.1", "MROak1.1", "MROak1.2", "MROak1…
## $ stand_ID <chr> "MRMOM", "MROak2", "MROak1", "MROak1", "MROak1", "MR…
## $ trap_number <chr> NA, "1", "1", "2", "3", "1", "2", "3", "1", "2", "3"…
## $ longitude <dbl> -73.59218, -73.58859, -73.58863, -73.58738, -73.5892…
## $ total_moth_count <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ complete <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ clean_complete <dbl> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ censored <int> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1…
## $ reason_incomplete <chr> NA, NA, NA, NA, NA, NA, NA, "Missing in field", NA, …
## $ X <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …stand_ID_filtered_1 <- stand_category_filtered %>%
separate(trap_name, into = c("stand_ID", "trap_number"), remove = FALSE, sep = "\\." ) %>%
glimpse()## Rows: 186
## Columns: 22
## $ patch_name <chr> "Mont_Royal", "Mont_Royal", "Mont_Royal", "Mont_Roya…
## $ Years <int> 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024, 2024…
## $ Year <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1…
## $ stand_type <chr> "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oak", "Oa…
## $ landscape_type <chr> "urban", "urban", "urban", "urban", "urban", "urban"…
## $ stand_category <chr> "E", "A", "C", "C", "C", "E", "E", "E", "E", "E", "E…
## $ Percent_Oak <dbl> 0.4, 0.8, 0.6, 0.6, 0.6, 0.4, 0.4, 0.4, 0.4, 0.4, 0.…
## $ Percent_Pine <dbl> 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.1, 0.1, 0.…
## $ Percent_Maple <dbl> 0.3, 0.1, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.0, 0.0, 0.…
## $ Forest_Type <chr> "oak_maple", "oak", "oak", "oak", "oak", "oak_maple"…
## $ stand_area_ha <dbl> 14.0, 6.4, 8.6, 8.6, 8.6, 11.7, 11.7, 11.7, 8.9, 8.9…
## $ forest_area_km2 <dbl> 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89…
## $ trap_name <chr> "MRMOM", "MROak2.1", "MROak1.1", "MROak1.2", "MROak1…
## $ stand_ID <chr> "MRMOM", "MROak2", "MROak1", "MROak1", "MROak1", "MR…
## $ trap_number <chr> NA, "1", "1", "2", "3", "1", "2", "3", "1", "2", "3"…
## $ longitude <dbl> -73.59218, -73.58859, -73.58863, -73.58738, -73.5892…
## $ total_moth_count <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ complete <int> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ clean_complete <dbl> 93, 15, 7, 20, 6, 96, 52, NA, 60, 20, 15, 17, 25, 3,…
## $ censored <int> 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1…
## $ reason_incomplete <chr> NA, NA, NA, NA, NA, NA, NA, "Missing in field", NA, …
## $ X <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …total_moths_2023_2024 <- sum(complete_2023_2024$clean_complete, na.rm = TRUE)
print(total_moths_2023_2024)## [1] 12456n_distinct(complete_2023_2024$trap_name)## [1] 255n_distinct(stand_ID_filtered$stand_ID)## [1] 155n_distinct(stand_ID_filtered$patch_name)## [1] 15#check to see the distribution of moth count data
hist(complete_2023_2024$clean_complete,
main = " ",
xlab = "Spongy moth count/trap",
ylab = "Frequency",
col = "darkblue",
border = "black")
hist(complete_2023_2024$clean_complete,
main = " ",
xlab = "Spongy moth count/trap",
ylab = "Frequency",
col = "blue",
border = "black",
breaks = 50) # increase this number to make bins smaller
#Calculate the mean and standard deviation of moth counts for each
#level of stand_category to see if there are differences.
# Checking unique combinations
table(stand_category_filtered$stand_category,
stand_category_filtered$patch_name)##
## Brownsburg Hatley Kenauk Mont_Gauvin/Glen Mont_Royal Mont_Saint_Bruno
## A 0 0 1 0 1 1
## B 0 0 0 0 0 0
## C 2 2 1 2 3 3
## D 0 0 0 0 0 0
## E 0 0 0 0 8 0
## F 0 0 0 0 0 0
## G 0 0 0 0 0 0
## H 0 0 0 0 0 0
## I 0 0 0 0 0 0
## J 0 0 0 0 0 0
## K 0 0 0 0 0 0
## L 2 0 0 0 1 0
## M 0 0 0 0 0 0
##
## Mont_Saint_Hilaire Montebello Notre_Dame_de_Bonsecours Oka Orford
## A 0 0 0 0 3
## B 0 2 0 0 0
## C 8 4 6 0 2
## D 3 0 0 0 0
## E 0 0 0 8 6
## F 2 0 5 0 5
## G 0 2 0 1 0
## H 4 0 1 8 0
## I 2 5 0 18 0
## J 0 0 0 0 4
## K 0 2 0 5 0
## L 1 0 0 2 7
## M 0 2 2 0 2
##
## Papineauville Parc_Michel_Chartrand Rigaud
## A 4 0 0
## B 0 0 0
## C 0 2 0
## D 0 0 7
## E 2 0 0
## F 0 0 5
## G 0 0 0
## H 2 0 3
## I 0 0 4
## J 0 0 3
## K 0 0 3
## L 0 0 2
## M 0 0 0table(stand_type_filtered$stand_type,
stand_type_filtered$patch_name)##
## Brownsburg Hatley Kenauk Mont_Gauvin/Glen Mont_Royal
## Oak 2 2 1 2 8
## Oak/Other 2 2 2 2 0
## Oak/Pine 0 0 1 0 4
## Other 0 2 2 2 0
## Pine 2 0 0 0 1
## Pine/Oak 0 0 0 0 0
##
## Mont_Saint_Bruno Mont_Saint_Hilaire Montebello
## Oak 3 8 6
## Oak/Other 5 3 2
## Oak/Pine 0 8 2
## Other 4 4 2
## Pine 0 1 5
## Pine/Oak 0 2 4
##
## Notre_Dame_de_Bonsecours Oka Orford Papineauville
## Oak 6 6 14 6
## Oak/Other 1 2 4 2
## Oak/Pine 5 9 2 2
## Other 2 0 4 2
## Pine 2 7 13 0
## Pine/Oak 1 18 0 0
##
## Parc_Michel_Chartrand Parc_Pointe_aux_Prairies Rigaud
## Oak 2 0 7
## Oak/Other 2 2 2
## Oak/Pine 0 0 8
## Other 2 4 0
## Pine 0 0 8
## Pine/Oak 0 0 4# Set the levels of 'stand_type' to ensure the correct order
stand_type_filtered$stand_type <- factor(stand_type_filtered$stand_type,
levels = c("Oak", "Oak/Pine", "Oak/Other",
"Pine/Oak", "Pine", "Other"))
# Create the 'stand-type by patch' table
table(stand_type_filtered$stand_type, stand_type_filtered$patch_name)##
## Brownsburg Hatley Kenauk Mont_Gauvin/Glen Mont_Royal
## Oak 2 2 1 2 8
## Oak/Pine 0 0 1 0 4
## Oak/Other 2 2 2 2 0
## Pine/Oak 0 0 0 0 0
## Pine 2 0 0 0 1
## Other 0 2 2 2 0
##
## Mont_Saint_Bruno Mont_Saint_Hilaire Montebello
## Oak 3 8 6
## Oak/Pine 0 8 2
## Oak/Other 5 3 2
## Pine/Oak 0 2 4
## Pine 0 1 5
## Other 4 4 2
##
## Notre_Dame_de_Bonsecours Oka Orford Papineauville
## Oak 6 6 14 6
## Oak/Pine 5 9 2 2
## Oak/Other 1 2 4 2
## Pine/Oak 1 18 0 0
## Pine 2 7 13 0
## Other 2 0 4 2
##
## Parc_Michel_Chartrand Parc_Pointe_aux_Prairies Rigaud
## Oak 2 0 7
## Oak/Pine 0 0 8
## Oak/Other 2 2 2
## Pine/Oak 0 0 4
## Pine 0 0 8
## Other 2 4 0# Create the contingency table
contingency_table <- table(stand_type_filtered$stand_type, stand_type_filtered$patch_name)
# Convert the table to a data frame
contingency_df <- as.data.frame(contingency_table)
# Heatmap stands by patch -------------------------------------------------
# Create a plot (heatmap) of the contingency table
contingency_df ## Var1 Var2 Freq
## 1 Oak Brownsburg 2
## 2 Oak/Pine Brownsburg 0
## 3 Oak/Other Brownsburg 2
## 4 Pine/Oak Brownsburg 0
## 5 Pine Brownsburg 2
## 6 Other Brownsburg 0
## 7 Oak Hatley 2
## 8 Oak/Pine Hatley 0
## 9 Oak/Other Hatley 2
## 10 Pine/Oak Hatley 0
## 11 Pine Hatley 0
## 12 Other Hatley 2
## 13 Oak Kenauk 1
## 14 Oak/Pine Kenauk 1
## 15 Oak/Other Kenauk 2
## 16 Pine/Oak Kenauk 0
## 17 Pine Kenauk 0
## 18 Other Kenauk 2
## 19 Oak Mont_Gauvin/Glen 2
## 20 Oak/Pine Mont_Gauvin/Glen 0
## 21 Oak/Other Mont_Gauvin/Glen 2
## 22 Pine/Oak Mont_Gauvin/Glen 0
## 23 Pine Mont_Gauvin/Glen 0
## 24 Other Mont_Gauvin/Glen 2
## 25 Oak Mont_Royal 8
## 26 Oak/Pine Mont_Royal 4
## 27 Oak/Other Mont_Royal 0
## 28 Pine/Oak Mont_Royal 0
## 29 Pine Mont_Royal 1
## 30 Other Mont_Royal 0
## 31 Oak Mont_Saint_Bruno 3
## 32 Oak/Pine Mont_Saint_Bruno 0
## 33 Oak/Other Mont_Saint_Bruno 5
## 34 Pine/Oak Mont_Saint_Bruno 0
## 35 Pine Mont_Saint_Bruno 0
## 36 Other Mont_Saint_Bruno 4
## 37 Oak Mont_Saint_Hilaire 8
## 38 Oak/Pine Mont_Saint_Hilaire 8
## 39 Oak/Other Mont_Saint_Hilaire 3
## 40 Pine/Oak Mont_Saint_Hilaire 2
## 41 Pine Mont_Saint_Hilaire 1
## 42 Other Mont_Saint_Hilaire 4
## 43 Oak Montebello 6
## 44 Oak/Pine Montebello 2
## 45 Oak/Other Montebello 2
## 46 Pine/Oak Montebello 4
## 47 Pine Montebello 5
## 48 Other Montebello 2
## 49 Oak Notre_Dame_de_Bonsecours 6
## 50 Oak/Pine Notre_Dame_de_Bonsecours 5
## 51 Oak/Other Notre_Dame_de_Bonsecours 1
## 52 Pine/Oak Notre_Dame_de_Bonsecours 1
## 53 Pine Notre_Dame_de_Bonsecours 2
## 54 Other Notre_Dame_de_Bonsecours 2
## 55 Oak Oka 6
## 56 Oak/Pine Oka 9
## 57 Oak/Other Oka 2
## 58 Pine/Oak Oka 18
## 59 Pine Oka 7
## 60 Other Oka 0
## 61 Oak Orford 14
## 62 Oak/Pine Orford 2
## 63 Oak/Other Orford 4
## 64 Pine/Oak Orford 0
## 65 Pine Orford 13
## 66 Other Orford 4
## 67 Oak Papineauville 6
## 68 Oak/Pine Papineauville 2
## 69 Oak/Other Papineauville 2
## 70 Pine/Oak Papineauville 0
## 71 Pine Papineauville 0
## 72 Other Papineauville 2
## 73 Oak Parc_Michel_Chartrand 2
## 74 Oak/Pine Parc_Michel_Chartrand 0
## 75 Oak/Other Parc_Michel_Chartrand 2
## 76 Pine/Oak Parc_Michel_Chartrand 0
## 77 Pine Parc_Michel_Chartrand 0
## 78 Other Parc_Michel_Chartrand 2
## 79 Oak Parc_Pointe_aux_Prairies 0
## 80 Oak/Pine Parc_Pointe_aux_Prairies 0
## 81 Oak/Other Parc_Pointe_aux_Prairies 2
## 82 Pine/Oak Parc_Pointe_aux_Prairies 0
## 83 Pine Parc_Pointe_aux_Prairies 0
## 84 Other Parc_Pointe_aux_Prairies 4
## 85 Oak Rigaud 7
## 86 Oak/Pine Rigaud 8
## 87 Oak/Other Rigaud 2
## 88 Pine/Oak Rigaud 4
## 89 Pine Rigaud 8
## 90 Other Rigaud 0ggplot(contingency_df, aes(x = Var1, y = Var2, fill = Freq)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "darkred") +
labs(x = "Stand Type", y = "Patch Name", fill = "Frequency") +
theme_minimal()
# Save the plot as an image (e.g., PNG)
#ggsave("contingency_table_heatmap.png", width = 7.5, height = 6)
#convert the data in the heatmap plot to a table format
wide_table <- contingency_df %>%
pivot_wider(names_from = Var1, values_from = Freq, values_fill = 0)
wide_table <- wide_table %>%
mutate(
Var2 = trimws(as.character(Var2)),
Var2 = dplyr::recode(Var2,
"Mont_Gauvin/Glen" = "Mont Gauvin",
"Mont_Royal" = "Mont Royal",
"Mont_Saint_Bruno" = "Mont Saint Bruno",
"Mont_Saint_Hilaire" = "Mont Saint Hilaire",
"Notre_Dame_de_Bonsecours" = "Notre Dame de Bonsecours",
"Oka" = "Parc Oka",
"Orford" = "Mont Orford",
"Parc_Michel_Chartrand" = "Parc Michel Chartrand",
"Parc_Pointe_aux_Prairies" = "Parc Pointe aux Prairies",
"Rigaud" = "Mont Rigaud",
))
##print(wide_table)
patch_order <- c("Papineauville", "Montebello", "Notre Dame de Bonsecours",
"Kenauk", "Brownsburg", "Mont Rigaud", "Parc Oka",
"Mont Royal",
"Parc Pointe aux Prairies", "Parc Michel Chartrand",
"Mont Saint Bruno", "Mont Saint Hilaire","Mont Gauvin",
"Mont Orford", "Hatley")
wide_table <- wide_table %>%
mutate(Var2 = factor(Var2, levels = patch_order)) %>%
arrange(Var2)
##print(wide_table)
gt_table <- wide_table %>%
gt() %>%
cols_label(
Var2 = "Patch Name"
) %>%
tab_options(
table.font.size = 12,
heading.title.font.size = 16,
heading.subtitle.font.size = 14
) %>%
cols_align(
align = "center",
columns = everything()
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels(everything())
) %>%
tab_style( # Custom striping on even-numbered rows
style = cell_fill(color = "#f9f9f0"), # You can change this color
locations = cells_body(rows = seq(2, nrow(wide_table), 2))
)
version$version.string## [1] "R version 4.4.3 (2025-02-28)"# Summary statistics for moth count by stand_category and stand_type
summary_stats <- stand_type_filtered %>%
group_by(stand_category, stand_type, patch_name) %>%
summarise(
mean_count = mean(total_moth_count, na.rm = TRUE),
sd_count = sd(total_moth_count, na.rm = TRUE),
var_count = var(total_moth_count, na.rm = TRUE),
count = n()
)
##print(summary_stats, n=22)
# Table to use -----------------------------------------------------------
# Summary statistics for moth count by stand_type, differentiating between
#traps set and traps with usable data
summary_stats_3 <- stand_type_filtered %>%
group_by(stand_type) %>%
summarise(
mean_count = mean(clean_complete, na.rm = TRUE),
sd_count = sd(clean_complete, na.rm = TRUE),
n_traps = n(),
n_obs = sum(!is.na (clean_complete))
)
print(summary_stats_3, n=22)## # A tibble: 6 × 5
## stand_type mean_count sd_count n_traps n_obs
## <fct> <dbl> <dbl> <int> <int>
## 1 Oak 77.0 69.8 73 57
## 2 Oak/Pine 64.4 52.6 41 30
## 3 Oak/Other 37.3 27.0 33 29
## 4 Pine/Oak 85.8 61.1 29 20
## 5 Pine 54.8 33.2 39 27
## 6 Other 36.9 26.8 30 27# Round numeric columns to 2 decimal places
summary_stats_3$mean_count <- round(summary_stats_3$mean_count, 2)
summary_stats_3$sd_count <- round(summary_stats_3$sd_count, 2)
summary_stats_3$count <- round(summary_stats_3$n_obs, 2)
# Summary statistics for moth count by patch
summary_stats_4 <- stand_type_filtered %>%
group_by(patch_name) %>%
summarise(
mean_count = mean(clean_complete, na.rm = TRUE),
sd_count = sd(clean_complete, na.rm = TRUE),
count = n()
)
print(summary_stats_4, n=22)## # A tibble: 15 × 4
## patch_name mean_count sd_count count
## <chr> <dbl> <dbl> <int>
## 1 "Brownsburg" 38.6 9.40 6
## 2 "Hatley" 74.8 18.4 6
## 3 "Kenauk" 26.8 20.4 6
## 4 "Mont_Gauvin/Glen " 34.1 11.6 6
## 5 "Mont_Royal" 35.5 32.0 13
## 6 "Mont_Saint_Bruno" 24.7 13.9 12
## 7 "Mont_Saint_Hilaire" 23.5 29.6 26
## 8 "Montebello" 54 44.8 21
## 9 "Notre_Dame_de_Bonsecours" 54.8 32.9 17
## 10 "Oka" 105. 63.6 42
## 11 "Orford" 88.0 67.9 37
## 12 "Papineauville" 71.0 56.2 12
## 13 "Parc_Michel_Chartrand" 70.2 17.4 6
## 14 "Parc_Pointe_aux_Prairies" 43.3 32.7 6
## 15 "Rigaud" 38.1 26.4 29# Summary statistics for moth count by stand_type and patch
# summary_stats_5 <- stand_type_filtered %>%
# Remove MOM row in 'complete' moth counts
moth_by_stand_summary_stats <- summary_stats_3[-c(1),]
print(moth_by_stand_summary_stats, n=22)## # A tibble: 5 × 6
## stand_type mean_count sd_count n_traps n_obs count
## <fct> <dbl> <dbl> <int> <int> <dbl>
## 1 Oak/Pine 64.4 52.6 41 30 30
## 2 Oak/Other 37.3 27.0 33 29 29
## 3 Pine/Oak 85.8 61.1 29 20 20
## 4 Pine 54.8 33.2 39 27 27
## 5 Other 36.9 26.8 30 27 27# Remove MOM row in 'clean_complete' moth counts
moth_by_stand_summary_stats_2 <- summary_stats_4[-c(1),]
print(moth_by_stand_summary_stats_2, n=22)## # A tibble: 14 × 4
## patch_name mean_count sd_count count
## <chr> <dbl> <dbl> <int>
## 1 "Hatley" 74.8 18.4 6
## 2 "Kenauk" 26.8 20.4 6
## 3 "Mont_Gauvin/Glen " 34.1 11.6 6
## 4 "Mont_Royal" 35.5 32.0 13
## 5 "Mont_Saint_Bruno" 24.7 13.9 12
## 6 "Mont_Saint_Hilaire" 23.5 29.6 26
## 7 "Montebello" 54 44.8 21
## 8 "Notre_Dame_de_Bonsecours" 54.8 32.9 17
## 9 "Oka" 105. 63.6 42
## 10 "Orford" 88.0 67.9 37
## 11 "Papineauville" 71.0 56.2 12
## 12 "Parc_Michel_Chartrand" 70.2 17.4 6
## 13 "Parc_Pointe_aux_Prairies" 43.3 32.7 6
## 14 "Rigaud" 38.1 26.4 29# Visualizations ----------------------------------------------------------
##separate Stand Type column so that we have a Stand ID for each stand in each patch
#stand_ID_filtered <- stand_type_filtered %>%
#Poisson, using all levels of data collection as a random effect, except Oak
#which is being fitted as a fixed effect
model_complete_poisson_oak <- glmer(
round(clean_complete) ~ (1|trap_name) +
(Percent_Oak) + (1|patch_name),
family =poisson(), data = stand_ID_filtered)
summary(model_complete_poisson_oak)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ (1 | trap_name) + (Percent_Oak) + (1 |
## patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1912.6 1925.6 -952.3 1904.6 186
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06666 -0.11714 0.00811 0.08998 0.36154
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.588 0.7668
## patch_name (Intercept) 0.262 0.5118
## Number of obs: 190, groups: trap_name, 190; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3637 0.1752 19.201 < 2e-16 ***
## Percent_Oak 0.7238 0.2630 2.752 0.00592 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Oak -0.525performance::check_overdispersion(model_complete_poisson_oak)## # Overdispersion test
##
## dispersion ratio = 0.142
## Pearson's Chi-Squared = 26.378
## p-value = 1performance::check_model(model_complete_poisson_oak)
#Heavily underdispersed, indicating that the moth counts are less variable
#than expected
#Run the same basic model, but have stands nested within patches as a
#random effects, first keeping in trap_name as a random effect and then
#removing it
model_complete_poisson_oak_nested <- glmer(
round(clean_complete) ~ Percent_Oak +
(1 | trap_name) +
(1 | patch_name/stand_ID), # nesting structure
family = poisson(),
data = stand_ID_filtered
)
summary(model_complete_poisson_oak_nested)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ Percent_Oak + (1 | trap_name) + (1 |
## patch_name/stand_ID)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1905.2 1921.5 -947.6 1895.2 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96329 -0.11714 0.00997 0.09144 0.37052
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.3346 0.5785
## stand_ID:patch_name (Intercept) 0.2975 0.5454
## patch_name (Intercept) 0.2328 0.4824
## Number of obs: 190, groups:
## trap_name, 190; stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3662 0.1785 18.862 <2e-16 ***
## Percent_Oak 0.6520 0.2997 2.175 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Oak -0.578performance::check_overdispersion(model_complete_poisson_oak_nested)## # Overdispersion test
##
## dispersion ratio = 0.146
## Pearson's Chi-Squared = 27.003
## p-value = 1performance::check_model(model_complete_poisson_oak_nested)
model_complete_poisson_oak_nested_1 <- glmer(
round(clean_complete) ~ Percent_Oak +
#(1 | trap_name) +
(1 | patch_name/stand_ID), # nesting structure
family = poisson(),
data = stand_ID_filtered
)
summary(model_complete_poisson_oak_nested_1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ Percent_Oak + (1 | patch_name/stand_ID)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 2887.1 2900.1 -1439.6 2879.1 186
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.0329 -0.7169 -0.0228 0.1485 13.3173
##
## Random effects:
## Groups Name Variance Std.Dev.
## stand_ID:patch_name (Intercept) 0.6051 0.7779
## patch_name (Intercept) 0.2183 0.4672
## Number of obs: 190, groups: stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3726 0.1777 18.982 <2e-16 ***
## Percent_Oak 0.6415 0.3083 2.081 0.0374 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Oak -0.596performance::check_overdispersion(model_complete_poisson_oak_nested_1)## # Overdispersion test
##
## dispersion ratio = 6.476
## Pearson's Chi-Squared = 1204.591
## p-value = < 0.001performance::check_model(model_complete_poisson_oak_nested_1)
#Try 3 levels of random effects
model_complete_poisson_oak_nested_2 <- glmer(
round(clean_complete) ~ Percent_Oak +
(1 | patch_name/stand_ID/trap_name), # nesting structure
family = poisson(),
data = stand_ID_filtered
)
summary(model_complete_poisson_oak_nested_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Oak + (1 | patch_name/stand_ID/trap_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1905.2 1921.5 -947.6 1895.2 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96329 -0.11714 0.00997 0.09144 0.37052
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name:stand_ID:patch_name (Intercept) 0.3346 0.5785
## stand_ID:patch_name (Intercept) 0.2975 0.5454
## patch_name (Intercept) 0.2328 0.4824
## Number of obs: 190, groups:
## trap_name:stand_ID:patch_name, 190; stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.3662 0.1785 18.862 <2e-16 ***
## Percent_Oak 0.6520 0.2997 2.175 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Oak -0.578performance::check_overdispersion(model_complete_poisson_oak_nested_2)## # Overdispersion test
##
## dispersion ratio = 0.146
## Pearson's Chi-Squared = 27.003
## p-value = 1performance::check_model(model_complete_poisson_oak_nested_2)
#Poisson, using all levels of data collection as a random effect, except Pine
#which is being fitted as a fixed effect
model_complete_poisson_pine <- glmer(
round(clean_complete) ~ (1|trap_name) +
(Percent_Pine) + (1|patch_name),
family =poisson(), data = stand_ID_filtered)
summary(model_complete_poisson_pine)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ (1 | trap_name) + (Percent_Pine) + (1 |
## patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1919.8 1932.8 -955.9 1911.8 186
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.94420 -0.11518 0.00723 0.08474 0.29169
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.6211 0.7881
## patch_name (Intercept) 0.2417 0.4916
## Number of obs: 190, groups: trap_name, 190; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5990 0.1499 24.015 <2e-16 ***
## Percent_Pine 0.1201 0.2866 0.419 0.675
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Pin -0.245performance::check_overdispersion(model_complete_poisson_pine)## # Overdispersion test
##
## dispersion ratio = 0.137
## Pearson's Chi-Squared = 25.519
## p-value = 1performance::check_model(model_complete_poisson_pine)
#Again, heavily underdispersed, indicating that the moth counts are
#less variable than expected
##Run the same basic model, but have stands nested within patches as a
#random effects, first keeping in trap_name as a random effect and then
#removing it
model_complete_poisson_pine_nested <- glmer(
round(clean_complete) ~ (Percent_Pine) +
(1 | trap_name) +
(1 | patch_name/stand_ID), # nesting structure
family = poisson(),
data = stand_ID_filtered)
summary(model_complete_poisson_pine_nested)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ (Percent_Pine) + (1 | trap_name) + (1 |
## patch_name/stand_ID)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1908.4 1924.6 -949.2 1898.4 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.83251 -0.13372 0.00920 0.08245 0.39874
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.3234 0.5687
## stand_ID:patch_name (Intercept) 0.3452 0.5875
## patch_name (Intercept) 0.2095 0.4577
## Number of obs: 190, groups:
## trap_name, 190; stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5375 0.1478 23.936 <2e-16 ***
## Percent_Pine 0.4124 0.3502 1.177 0.239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Pin -0.291performance::check_overdispersion(model_complete_poisson_pine_nested)## # Overdispersion test
##
## dispersion ratio = 0.142
## Pearson's Chi-Squared = 26.303
## p-value = 1performance::check_model(model_complete_poisson_pine_nested)
model_complete_poisson_pine_nested_1 <- glmer(
round(clean_complete) ~ (Percent_Pine) +
#(1 | trap_name) +
(1 | patch_name/stand_ID), # nesting structure
family = poisson(),
data = stand_ID_filtered)
summary(model_complete_poisson_pine_nested_1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ (Percent_Pine) + (1 | patch_name/stand_ID)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 2889.1 2902.1 -1440.6 2881.1 186
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.0307 -0.6664 -0.0092 0.1470 13.3232
##
## Random effects:
## Groups Name Variance Std.Dev.
## stand_ID:patch_name (Intercept) 0.6281 0.7925
## patch_name (Intercept) 0.1994 0.4466
## Number of obs: 190, groups: stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5296 0.1451 24.331 <2e-16 ***
## Percent_Pine 0.5442 0.3650 1.491 0.136
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Pin -0.287performance::check_overdispersion(model_complete_poisson_pine_nested_1)## # Overdispersion test
##
## dispersion ratio = 6.474
## Pearson's Chi-Squared = 1204.143
## p-value = < 0.001performance::check_model(model_complete_poisson_pine_nested_1)
model_complete_poisson_pine_nested_2 <- glmer(
round(clean_complete) ~ (Percent_Pine) +
(1 | patch_name/stand_ID/trap_name), # nesting structure
family = poisson(),
data = stand_ID_filtered)
summary(model_complete_poisson_pine_nested_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ (Percent_Pine) + (1 | patch_name/stand_ID/trap_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1908.4 1924.6 -949.2 1898.4 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.83251 -0.13372 0.00920 0.08245 0.39874
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name:stand_ID:patch_name (Intercept) 0.3234 0.5687
## stand_ID:patch_name (Intercept) 0.3452 0.5875
## patch_name (Intercept) 0.2095 0.4577
## Number of obs: 190, groups:
## trap_name:stand_ID:patch_name, 190; stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.5375 0.1478 23.936 <2e-16 ***
## Percent_Pine 0.4124 0.3502 1.177 0.239
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Pin -0.291performance::check_overdispersion(model_complete_poisson_pine_nested_2)## # Overdispersion test
##
## dispersion ratio = 0.142
## Pearson's Chi-Squared = 26.303
## p-value = 1performance::check_model(model_complete_poisson_pine_nested_2)
#model for Oak and Pine together
model_both <- glmer(round(clean_complete) ~ Percent_Pine + Percent_Oak +
(1 | trap_name) + (1 | patch_name),
data = stand_ID_filtered, family = poisson)
summary(model_both)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | trap_name) +
## (1 | patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1909.4 1925.6 -949.7 1899.4 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06551 -0.12913 0.00495 0.09740 0.32161
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5761 0.7590
## patch_name (Intercept) 0.2482 0.4981
## Number of obs: 190, groups: trap_name, 190; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1338 0.1995 15.705 < 2e-16 ***
## Percent_Pine 0.7514 0.3292 2.283 0.022448 *
## Percent_Oak 1.1112 0.3113 3.569 0.000358 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_P
## Percent_Pin -0.509
## Percent_Oak -0.660 0.546performance::check_overdispersion(model_both)## # Overdispersion test
##
## dispersion ratio = 0.141
## Pearson's Chi-Squared = 26.063
## p-value = 1#Again, heavily underdispersed, indicating that the moth counts are
#less variable than expected
##Run the same basic model, but have stands nested within patches as a
#random effects, first keeping in trap_name as a random effect and then
#removing it
model_both_nested <- glmer(
round(clean_complete) ~ Percent_Pine + Percent_Oak +
(1 | trap_name) +
(1 | patch_name/stand_ID), # nesting structure
family = poisson(),
data = stand_ID_filtered)
summary(model_both_nested)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | trap_name) +
## (1 | patch_name/stand_ID)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1901.0 1920.5 -944.5 1889.0 184
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96859 -0.13088 0.00960 0.09818 0.38032
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.3268 0.5716
## stand_ID:patch_name (Intercept) 0.2859 0.5347
## patch_name (Intercept) 0.2239 0.4732
## Number of obs: 190, groups:
## trap_name, 190; stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1157 0.2023 15.399 < 2e-16 ***
## Percent_Pine 0.9440 0.3774 2.501 0.01238 *
## Percent_Oak 1.0561 0.3368 3.136 0.00171 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_P
## Percent_Pin -0.498
## Percent_Oak -0.678 0.473performance::check_overdispersion(model_both_nested)## # Overdispersion test
##
## dispersion ratio = 0.146
## Pearson's Chi-Squared = 26.852
## p-value = 1model_both_nested_1 <- glmer(
round(clean_complete) ~ Percent_Pine + Percent_Oak +
#(1 | trap_name) +
(1 | patch_name/stand_ID), # nesting structure
family = poisson(),
data = stand_ID_filtered)
summary(model_both_nested_1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | patch_name/stand_ID)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 2881.9 2898.1 -1435.9 2871.9 185
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -9.0299 -0.7789 -0.0275 0.1466 13.3255
##
## Random effects:
## Groups Name Variance Std.Dev.
## stand_ID:patch_name (Intercept) 0.5731 0.7570
## patch_name (Intercept) 0.2117 0.4601
## Number of obs: 190, groups: stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1210 0.1980 15.761 < 2e-16 ***
## Percent_Pine 1.0561 0.3896 2.711 0.00671 **
## Percent_Oak 1.0370 0.3349 3.097 0.00196 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_P
## Percent_Pin -0.476
## Percent_Oak -0.677 0.438performance::check_overdispersion(model_both_nested_1)## # Overdispersion test
##
## dispersion ratio = 6.509
## Pearson's Chi-Squared = 1204.198
## p-value = < 0.001model_both_nested_2 <- glmer(
round(clean_complete) ~ Percent_Pine + Percent_Oak +
(1 | patch_name/stand_ID/trap_name), # nesting structure
family = poisson(),
data = stand_ID_filtered)
summary(model_both_nested_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | patch_name/stand_ID/trap_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1901.0 1920.5 -944.5 1889.0 184
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96859 -0.13088 0.00960 0.09818 0.38032
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name:stand_ID:patch_name (Intercept) 0.3268 0.5716
## stand_ID:patch_name (Intercept) 0.2859 0.5347
## patch_name (Intercept) 0.2239 0.4732
## Number of obs: 190, groups:
## trap_name:stand_ID:patch_name, 190; stand_ID:patch_name, 131; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1157 0.2023 15.399 < 2e-16 ***
## Percent_Pine 0.9440 0.3774 2.501 0.01238 *
## Percent_Oak 1.0561 0.3368 3.136 0.00171 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_P
## Percent_Pin -0.498
## Percent_Oak -0.678 0.473performance::check_overdispersion(model_both_nested_2)## # Overdispersion test
##
## dispersion ratio = 0.146
## Pearson's Chi-Squared = 26.852
## p-value = 1#model for Oak and Pine together, adding an interaction of oak & pine
model_both_2 <- glmer(round(clean_complete) ~ Percent_Pine + Percent_Oak +
(1 | trap_name) + (1 | patch_name) +
(Percent_Pine * Percent_Oak),
data = stand_ID_filtered, family = poisson)
summary(model_both_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Pine + Percent_Oak + (1 | trap_name) +
## (1 | patch_name) + (Percent_Pine * Percent_Oak)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1908.7 1928.2 -948.4 1896.7 184
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.04594 -0.13354 0.01572 0.09499 0.35363
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5692 0.7544
## patch_name (Intercept) 0.2484 0.4984
## Number of obs: 190, groups: trap_name, 190; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.1521 0.1993 15.818 < 2e-16 ***
## Percent_Pine 0.4895 0.3640 1.345 0.17863
## Percent_Oak 0.9508 0.3249 2.926 0.00343 **
## Percent_Pine:Percent_Oak 2.7525 1.6875 1.631 0.10286
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_P Prcn_O
## Percent_Pin -0.479
## Percent_Oak -0.643 0.601
## Prcnt_P:P_O 0.051 -0.437 -0.299performance::check_overdispersion(model_both_2)## # Overdispersion test
##
## dispersion ratio = 0.140
## Pearson's Chi-Squared = 25.805
## p-value = 1ggplot(stand_ID_filtered, aes(x = Percent_Oak, y = clean_complete)) +
geom_point(alpha = 0.5) +
stat_smooth(method = "glm", method.args = list(family = "poisson"),
se = TRUE, color = "darkgreen") +
labs(
title = "Effect of Percent Oak on Moth Counts",
x = "% Oak Cover",
y = "Moth Count"
) +
theme_minimal()
# Checking Maple ----------------------------------------------------------
#Poisson, using all levels of data collection as a random effect, except Maple
#which is being fitted as a fixed effect
model_complete_poisson_maple <- glmer(
round(clean_complete) ~ (1|trap_name) +
(Percent_Maple) + (1|patch_name),
family =poisson(), data = stand_ID_filtered)
summary(model_complete_poisson_maple)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ (1 | trap_name) + (Percent_Maple) + (1 |
## patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1872.4 1885.3 -932.2 1864.4 182
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.92817 -0.13025 0.00652 0.09274 0.30691
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5817 0.7627
## patch_name (Intercept) 0.2590 0.5090
## Number of obs: 186, groups: trap_name, 186; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.7646 0.1869 20.144 <2e-16 ***
## Percent_Maple -0.4860 0.4528 -1.073 0.283
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## Percent_Mpl -0.604performance::check_overdispersion(model_complete_poisson_maple)## # Overdispersion test
##
## dispersion ratio = 0.127
## Pearson's Chi-Squared = 23.172
## p-value = 1performance::check_model(model_complete_poisson_maple)
#model for Oak and Maple together
model_oak_maple <- glmer(round(clean_complete) ~ Percent_Maple + Percent_Oak +
(1 | trap_name) + (1 | patch_name),
data = stand_ID_filtered, family = poisson)
summary(model_oak_maple)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Maple + Percent_Oak + (1 | trap_name) +
## (1 | patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1863.7 1879.8 -926.8 1853.7 181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05891 -0.12810 0.00800 0.08859 0.38383
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5399 0.7348
## patch_name (Intercept) 0.2856 0.5344
## Number of obs: 186, groups: trap_name, 186; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4304 0.2140 16.028 < 2e-16 ***
## Percent_Maple -0.3332 0.4421 -0.754 0.451113
## Percent_Oak 0.8734 0.2622 3.331 0.000864 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_M
## Percent_Mpl -0.556
## Percent_Oak -0.470 0.098performance::check_overdispersion(model_oak_maple)## # Overdispersion test
##
## dispersion ratio = 0.132
## Pearson's Chi-Squared = 23.921
## p-value = 1#model for Oak and Maple together, adding an interaction of oak & maple
model_oak_maple_2 <- glmer(round(clean_complete) ~ Percent_Maple +
Percent_Oak +
(1 | trap_name) + (1 | patch_name) +
(Percent_Maple * Percent_Oak),
data = stand_ID_filtered, family = poisson)
summary(model_oak_maple_2)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ Percent_Maple + Percent_Oak + (1 | trap_name) +
## (1 | patch_name) + (Percent_Maple * Percent_Oak)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1865.6 1884.9 -926.8 1853.6 180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.06075 -0.12872 0.00638 0.08843 0.37974
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5401 0.7349
## patch_name (Intercept) 0.2842 0.5331
## Number of obs: 186, groups: trap_name, 186; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.4649 0.2400 14.440 <2e-16 ***
## Percent_Maple -0.4672 0.6119 -0.763 0.4452
## Percent_Oak 0.7598 0.4444 1.710 0.0873 .
## Percent_Maple:Percent_Oak 0.5046 1.5930 0.317 0.7514
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_M Prcn_O
## Percent_Mpl -0.672
## Percent_Oak -0.614 0.600
## Prcnt_M:P_O 0.454 -0.691 -0.807performance::check_overdispersion(model_oak_maple_2)## # Overdispersion test
##
## dispersion ratio = 0.133
## Pearson's Chi-Squared = 23.914
## p-value = 1#model for Oak, Maple, and Pine together
model_all <- glmer(round(clean_complete) ~ Percent_Maple + Percent_Oak +
Percent_Pine + (1 | trap_name) + (1 | patch_name),
data = stand_ID_filtered, family = poisson)
summary(model_all)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ Percent_Maple + Percent_Oak + Percent_Pine +
## (1 | trap_name) + (1 | patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1860.8 1880.2 -924.4 1848.8 180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.07728 -0.14342 0.00230 0.09453 0.37783
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5308 0.7286
## patch_name (Intercept) 0.2614 0.5112
## Number of obs: 186, groups: trap_name, 186; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.9191 0.3126 9.338 < 2e-16 ***
## Percent_Maple 0.4709 0.5687 0.828 0.4076
## Percent_Oak 1.4067 0.3558 3.954 7.69e-05 ***
## Percent_Pine 0.9371 0.4249 2.206 0.0274 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_M Prcn_O
## Percent_Mpl -0.765
## Percent_Oak -0.741 0.491
## Percent_Pin -0.744 0.639 0.682performance::check_overdispersion(model_all)## # Overdispersion test
##
## dispersion ratio = 0.131
## Pearson's Chi-Squared = 23.652
## p-value = 1## All Variables -----------------------------------------------------------
##Fitting a poisson model with all of the possible response variables, to
#explore a correlation between all possible variables and moth counts
all_variables_poisson <- glmer(round(clean_complete) ~ (1|trap_name) +
Percent_Oak + Percent_Pine + landscape_type + longitude +
forest_area_km2 + stand_area_ha + (1|patch_name),
family = poisson(), data = stand_ID_filtered)
summary(all_variables_poisson)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ (1 | trap_name) + Percent_Oak + Percent_Pine +
## landscape_type + longitude + forest_area_km2 + stand_area_ha +
## (1 | patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1909.0 1941.5 -944.5 1889.0 180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87569 -0.13572 0.00032 0.09462 0.35788
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5393 0.7344
## patch_name (Intercept) 0.3123 0.5588
## Number of obs: 190, groups: trap_name, 190; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.6680214 14.0675971 0.474 0.63550
## Percent_Oak 0.9611935 0.3182652 3.020 0.00253 **
## Percent_Pine 0.5048426 0.3447402 1.464 0.14308
## landscape_typeforest 0.9008992 0.3873226 2.326 0.02002 *
## landscape_typeurban -0.2759427 0.3971989 -0.695 0.48723
## longitude 0.0466809 0.1914203 0.244 0.80733
## forest_area_km2 -0.0005727 0.0003957 -1.447 0.14777
## stand_area_ha -0.0134213 0.0193206 -0.695 0.48727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_O Prcn_P lndscp_typf lndscp_typr longtd frs__2
## Percent_Oak 0.181
## Percent_Pin 0.159 0.597
## lndscp_typf -0.236 -0.168 -0.220
## lndscp_typr 0.033 0.139 0.182 0.260
## longitude 1.000 0.190 0.166 -0.230 0.044
## forst_r_km2 0.493 0.123 0.134 -0.499 0.070 0.497
## stand_are_h -0.041 -0.048 -0.018 -0.028 0.024 -0.030 0.075
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.243731 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?AIC(all_variables_poisson)## [1] 1908.988print(summary(all_variables_poisson))## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## round(clean_complete) ~ (1 | trap_name) + Percent_Oak + Percent_Pine +
## landscape_type + longitude + forest_area_km2 + stand_area_ha +
## (1 | patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1909.0 1941.5 -944.5 1889.0 180
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.87569 -0.13572 0.00032 0.09462 0.35788
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5393 0.7344
## patch_name (Intercept) 0.3123 0.5588
## Number of obs: 190, groups: trap_name, 190; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.6680214 14.0675971 0.474 0.63550
## Percent_Oak 0.9611935 0.3182652 3.020 0.00253 **
## Percent_Pine 0.5048426 0.3447402 1.464 0.14308
## landscape_typeforest 0.9008992 0.3873226 2.326 0.02002 *
## landscape_typeurban -0.2759427 0.3971989 -0.695 0.48723
## longitude 0.0466809 0.1914203 0.244 0.80733
## forest_area_km2 -0.0005727 0.0003957 -1.447 0.14777
## stand_area_ha -0.0134213 0.0193206 -0.695 0.48727
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_O Prcn_P lndscp_typf lndscp_typr longtd frs__2
## Percent_Oak 0.181
## Percent_Pin 0.159 0.597
## lndscp_typf -0.236 -0.168 -0.220
## lndscp_typr 0.033 0.139 0.182 0.260
## longitude 1.000 0.190 0.166 -0.230 0.044
## forst_r_km2 0.493 0.123 0.134 -0.499 0.070 0.497
## stand_are_h -0.041 -0.048 -0.018 -0.028 0.024 -0.030 0.075
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.243731 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?tab_model(all_variables_poisson,
show.ci = FALSE, # hide CI since you only want SE & p
show.stat = TRUE, # show test statistics
p.style = "numeric") # show numeric p-values| round(clean complete) | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | Statistic | p |
| (Intercept) | 786.84 | 0.47 | 0.636 |
| Percent Oak | 2.61 | 3.02 | 0.003 |
| Percent Pine | 1.66 | 1.46 | 0.143 |
| landscape typeforest | 2.46 | 2.33 | 0.020 |
| landscape type [urban] | 0.76 | -0.69 | 0.487 |
| longitude | 1.05 | 0.24 | 0.807 |
| forest area km2 | 1.00 | -1.45 | 0.148 |
| stand area ha | 0.99 | -0.69 | 0.487 |
| Random Effects | |||
| σ2 | 0.56 | ||
| τ00 trap_name | 0.54 | ||
| τ00 patch_name | 0.31 | ||
| ICC | 0.36 | ||
| N trap_name | 190 | ||
| N patch_name | 15 | ||
| Observations | 190 | ||
| Marginal R2 / Conditional R2 | 0.211 / 0.495 | ||
#checking for collinearity between variables
#using the 'performance' package, which is built for linear
#mixed models (glmer)
collinearity <- check_collinearity(all_variables_poisson)
print(collinearity)## # Check for Multicollinearity
##
## Low Correlation
##
## Term VIF VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
## Percent_Oak 1.58 [1.35, 1.96] 1.26 0.63 [0.51, 0.74]
## Percent_Pine 1.65 [1.40, 2.04] 1.28 0.61 [0.49, 0.71]
## landscape_type 1.57 [1.34, 1.94] 1.25 0.64 [0.51, 0.74]
## longitude 1.37 [1.20, 1.70] 1.17 0.73 [0.59, 0.84]
## forest_area_km2 1.78 [1.50, 2.22] 1.33 0.56 [0.45, 0.67]
## stand_area_ha 1.01 [1.00, 180.64] 1.01 0.99 [0.01, 1.00]# Convert to a regular data frame
col_df <- as.data.frame(collinearity)
print(col_df)## Term VIF VIF_CI_low VIF_CI_high SE_factor Tolerance
## 1 Percent_Oak 1.581963 1.353069 1.959249 1.257761 0.6321259
## 2 Percent_Pine 1.647135 1.401846 2.042150 1.283408 0.6071148
## 3 landscape_type 1.570228 1.344309 1.944384 1.253087 0.6368502
## 4 longitude 1.370329 1.196894 1.696532 1.170610 0.7297520
## 5 forest_area_km2 1.782123 1.503393 2.215186 1.334962 0.5611286
## 6 stand_area_ha 1.014701 1.000001 180.642302 1.007324 0.9855115
## Tolerance_CI_low Tolerance_CI_high
## 1 0.510399526 0.7390606
## 2 0.489680063 0.7133451
## 3 0.514301733 0.7438766
## 4 0.589437855 0.8354955
## 5 0.451429289 0.6651623
## 6 0.005535802 0.9999988# checking for co-variation between numeric predictors
#in a Pearson correlation matrix
numeric_vars <- stand_ID_filtered[, c("Percent_Oak", "Percent_Pine", "longitude",
"forest_area_km2", "stand_area_ha")]
# Correlation matrix
cor(numeric_vars, use = "complete.obs")## Percent_Oak Percent_Pine longitude forest_area_km2
## Percent_Oak 1.00000000 -0.52500780 -0.15613141 -0.01062395
## Percent_Pine -0.52500780 1.00000000 -0.11765413 0.05010221
## longitude -0.15613141 -0.11765413 1.00000000 -0.41050141
## forest_area_km2 -0.01062395 0.05010221 -0.41050141 1.00000000
## stand_area_ha -0.03911951 0.09618211 0.03658612 -0.24775945
## stand_area_ha
## Percent_Oak -0.03911951
## Percent_Pine 0.09618211
## longitude 0.03658612
## forest_area_km2 -0.24775945
## stand_area_ha 1.00000000library(gt)
cor_mat <- cor(numeric_vars, use = "complete.obs")
cor_df <- as.data.frame(round(cor_mat, 3))
cor_df |> gt::gt() |> gt::tab_caption("Correlation Matrix")| Percent_Oak | Percent_Pine | longitude | forest_area_km2 | stand_area_ha |
|---|---|---|---|---|
| 1.000 | -0.525 | -0.156 | -0.011 | -0.039 |
| -0.525 | 1.000 | -0.118 | 0.050 | 0.096 |
| -0.156 | -0.118 | 1.000 | -0.411 | 0.037 |
| -0.011 | 0.050 | -0.411 | 1.000 | -0.248 |
| -0.039 | 0.096 | 0.037 | -0.248 | 1.000 |
#write.csv(cor_df, "correlation_matrix.csv", row.names = TRUE)
#checking for interaction between oak and landscape type
interaction_model <- lm(clean_complete ~ Percent_Oak * landscape_type +
Percent_Pine + longitude + forest_area_km2 + stand_area_ha,
data = stand_ID_filtered)
summary(interaction_model)##
## Call:
## lm(formula = clean_complete ~ Percent_Oak * landscape_type +
## Percent_Pine + longitude + forest_area_km2 + stand_area_ha,
## data = stand_ID_filtered)
##
## Residuals:
## Min 1Q Median 3Q Max
## -68.68 -32.55 -12.11 22.52 287.54
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 281.52445 452.00476 0.623 0.53418
## Percent_Oak 43.51391 26.86459 1.620 0.10704
## landscape_typeforest 0.91995 17.39731 0.053 0.95789
## landscape_typeurban 15.38240 23.83683 0.645 0.51954
## Percent_Pine 56.61962 21.26213 2.663 0.00845 **
## longitude 3.42628 6.17707 0.555 0.57980
## forest_area_km2 -0.02668 0.01540 -1.732 0.08493 .
## stand_area_ha 0.56936 1.12480 0.506 0.61334
## Percent_Oak:landscape_typeforest 65.12563 35.49584 1.835 0.06819 .
## Percent_Oak:landscape_typeurban -66.97300 54.74304 -1.223 0.22278
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 51.2 on 180 degrees of freedom
## (55 observations deleted due to missingness)
## Multiple R-squared: 0.1451, Adjusted R-squared: 0.1024
## F-statistic: 3.396 on 9 and 180 DF, p-value: 0.0007096AIC(interaction_model)## [1] 2046.479#checking for interaction between oak and landscape type, with a glm model
interaction_model_glm <- glmer(round(clean_complete) ~
Percent_Oak * landscape_type +
Percent_Pine +
(1 | trap_name) + (1 | patch_name),
data = stand_ID_filtered, family = poisson)
summary(interaction_model_glm)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: round(clean_complete) ~ Percent_Oak * landscape_type + Percent_Pine +
## (1 | trap_name) + (1 | patch_name)
## Data: stand_ID_filtered
##
## AIC BIC logLik -2*log(L) df.resid
## 1898.4 1927.6 -940.2 1880.4 181
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.00781 -0.15339 0.00673 0.09710 0.37024
##
## Random effects:
## Groups Name Variance Std.Dev.
## trap_name (Intercept) 0.5103 0.7143
## patch_name (Intercept) 0.3576 0.5980
## Number of obs: 190, groups: trap_name, 190; patch_name, 15
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.9314 0.2906 10.089 < 2e-16 ***
## Percent_Oak 1.1848 0.4106 2.885 0.00391 **
## landscape_typeforest 0.5073 0.4099 1.238 0.21587
## landscape_typeurban 0.6694 0.5066 1.321 0.18642
## Percent_Pine 0.8037 0.3480 2.310 0.02090 *
## Percent_Oak:landscape_typeforest 0.4032 0.5333 0.756 0.44969
## Percent_Oak:landscape_typeurban -2.8278 0.9080 -3.114 0.00184 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Prcn_O lndscp_typf lndscp_typr Prcn_P
## Percent_Oak -0.535
## lndscp_typf -0.526 0.199
## lndscp_typr -0.556 0.328 0.279
## Percent_Pin -0.296 0.365 -0.263 0.230
## Prcnt_Ok:lndscp_typf 0.268 -0.594 -0.493 -0.125 0.201
## Prcnt_Ok:lndscp_typr 0.244 -0.447 -0.134 -0.573 -0.150
## Prcnt_Ok:lndscp_typf
## Percent_Oak
## lndscp_typf
## lndscp_typr
## Percent_Pin
## Prcnt_Ok:lndscp_typf
## Prcnt_Ok:lndscp_typr 0.267AIC(interaction_model_glm)## [1] 1898.357performance::check_overdispersion(interaction_model_glm)## # Overdispersion test
##
## dispersion ratio = 0.141
## Pearson's Chi-Squared = 25.499
## p-value = 1performance::check_model(interaction_model_glm)
tab_model(interaction_model_glm,
show.ci = FALSE, # hide CI since you only want SE & p
show.stat = TRUE, # show test statistics
p.style = "numeric") # show numeric p-values| round(clean complete) | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | Statistic | p |
| (Intercept) | 18.75 | 10.09 | <0.001 |
| Percent Oak | 3.27 | 2.89 | 0.004 |
| landscape typeforest | 1.66 | 1.24 | 0.216 |
| landscape type [urban] | 1.95 | 1.32 | 0.186 |
| Percent Pine | 2.23 | 2.31 | 0.021 |
|
Percent Oak × landscape typeforest |
1.50 | 0.76 | 0.450 |
|
Percent Oak × landscape type [urban] |
0.06 | -3.11 | 0.002 |
| Random Effects | |||
| σ2 | 0.53 | ||
| τ00 trap_name | 0.51 | ||
| τ00 patch_name | 0.36 | ||
| ICC | 0.40 | ||
| N trap_name | 190 | ||
| N patch_name | 15 | ||
| Observations | 190 | ||
| Marginal R2 / Conditional R2 | 0.204 / 0.526 | ||
ggplot(stand_ID_filtered, aes(x = Percent_Oak, y = clean_complete, color = landscape_type)) +
geom_point(alpha = 0.6) + # Add raw data points with transparency
geom_smooth(method = "lm", se = FALSE) + # Add regression lines by group
theme_minimal() +
labs(
title = " ",
x = "Percent Oak",
y = "Moth Counts",
color = "Landscape Type"
) +
scale_color_viridis_d(option = "D") # You can also try options "A", "B", "C", "E"
# Heatmap oak/pine spectrum -----------------------------------------------
# Create the 'oak/pine spectrum' table
table(stand_category_filtered$Percent_Oak,
stand_category_filtered$Percent_Pine)##
## 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
## 0 0 0 0 0 1 4 4 8 6
## 0.1 0 0 0 0 0 2 3 7 0
## 0.2 0 0 0 1 5 1 0 0 0
## 0.3 0 0 10 6 12 12 2 0 0
## 0.4 5 12 12 5 1 0 0 0 0
## 0.5 14 2 0 0 2 0 0 0 0
## 0.6 15 3 0 0 0 0 0 0 0
## 0.7 12 7 0 0 0 0 0 0 0
## 0.8 10 2 0 0 0 0 0 0 0# Create the contingency table
contingency_table <- table(stand_type_filtered$Percent_Oak,
stand_type_filtered$Percent_Pine)
# Convert the table to a data frame
contingency_df <- as.data.frame(contingency_table)
# Create a plot (heatmap) of the contingency table
contingency_df ## Var1 Var2 Freq
## 1 0 0 19
## 2 0.1 0 4
## 3 0.2 0 2
## 4 0.3 0 13
## 5 0.4 0 20
## 6 0.5 0 14
## 7 0.6 0 15
## 8 0.7 0 12
## 9 0.8 0 10
## 10 0 0.1 0
## 11 0.1 0.1 2
## 12 0.2 0.1 0
## 13 0.3 0.1 0
## 14 0.4 0.1 12
## 15 0.5 0.1 2
## 16 0.6 0.1 3
## 17 0.7 0.1 7
## 18 0.8 0.1 2
## 19 0 0.2 3
## 20 0.1 0.2 0
## 21 0.2 0.2 0
## 22 0.3 0.2 11
## 23 0.4 0.2 12
## 24 0.5 0.2 0
## 25 0.6 0.2 0
## 26 0.7 0.2 0
## 27 0.8 0.2 0
## 28 0 0.3 0
## 29 0.1 0.3 0
## 30 0.2 0.3 1
## 31 0.3 0.3 6
## 32 0.4 0.3 5
## 33 0.5 0.3 0
## 34 0.6 0.3 0
## 35 0.7 0.3 0
## 36 0.8 0.3 0
## 37 0 0.4 1
## 38 0.1 0.4 0
## 39 0.2 0.4 5
## 40 0.3 0.4 12
## 41 0.4 0.4 1
## 42 0.5 0.4 2
## 43 0.6 0.4 0
## 44 0.7 0.4 0
## 45 0.8 0.4 0
## 46 0 0.5 4
## 47 0.1 0.5 2
## 48 0.2 0.5 1
## 49 0.3 0.5 12
## 50 0.4 0.5 0
## 51 0.5 0.5 0
## 52 0.6 0.5 0
## 53 0.7 0.5 0
## 54 0.8 0.5 0
## 55 0 0.6 4
## 56 0.1 0.6 3
## 57 0.2 0.6 0
## 58 0.3 0.6 2
## 59 0.4 0.6 0
## 60 0.5 0.6 0
## 61 0.6 0.6 0
## 62 0.7 0.6 0
## 63 0.8 0.6 0
## 64 0 0.7 8
## 65 0.1 0.7 7
## 66 0.2 0.7 0
## 67 0.3 0.7 0
## 68 0.4 0.7 0
## 69 0.5 0.7 0
## 70 0.6 0.7 0
## 71 0.7 0.7 0
## 72 0.8 0.7 0
## 73 0 0.8 6
## 74 0.1 0.8 0
## 75 0.2 0.8 0
## 76 0.3 0.8 0
## 77 0.4 0.8 0
## 78 0.5 0.8 0
## 79 0.6 0.8 0
## 80 0.7 0.8 0
## 81 0.8 0.8 0ggplot(contingency_df, aes(x = Var1, y = Var2, fill = Freq)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "darkred") +
labs(x = "Proportion Oak", y = "Proportion Pine", fill = "Frequency") +
theme_minimal()
# Save the plot as an image (e.g., PNG)
#ggsave("oak_pine_heatmap.png", width = 7.5, height = 6)
# Citations ---------------------------------------------------------------
## @Manual{,
## title = {R: A Language and Environment for Statistical Computing},
## author = {{R Core Team}},
## organization = {R Foundation for Statistical Computing},
## address = {Vienna, Austria},
## year = {2021},
## url = {https://www.R-project.org/},
## }
citation("lme4")## To cite lme4 in publications use:
##
## Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015).
## Fitting Linear Mixed-Effects Models Using lme4. Journal of
## Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {Fitting Linear Mixed-Effects Models Using {lme4}},
## author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
## journal = {Journal of Statistical Software},
## year = {2015},
## volume = {67},
## number = {1},
## pages = {1--48},
## doi = {10.18637/jss.v067.i01},
## }xfun::Rscript_call( rmarkdown::render, list( input = “Report/2023_2024_All_Data_Oct_25.Rmd”, output_format = “html_document” ), fail = TRUE )
xfun::Rscript_call( rmarkdown::render, list(input = “Report/2023_2024_All_Data_Oct_25.Rmd”, output_format = “html_document”) )
xfun::Rscript_call( rmarkdown::render, list(input = “Report/2023_2024_All_Data_Oct_25.Rmd”, output_format = “pdf_document”) )